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The present study investigates the spatio-temporal variabiUty in the dynamics of self-sustained 
supersonic reaction waves propagating through an excitable medium. The model is an extension 
of Fickett's detonation model with a state dependent energy addition term. Stable and pulsating 
supersonic waves are predicted. With increasing sensitivity of the reaction rate, the reaction wave 
transits from steady propagation to stable limit cycles and eventually to chaos through the classical 
Feigenbaum route. The physical pulsation mechanism is explained by the coherence between internal 
wave motion and energy release. The results obtained clarify the physical origin of detonation wave 
instability in chemical detonations previously observed experimentally. 
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Self-sustained waves propagating through excitable 
media can give rise to complex patterns (e.g., spirals, 
cells, target shapes) and generally show large spatio- 
temporal variability in their dynamics. While sub- 
sonic wave propagation is usually modelled by reaction- 
diffusion coupling, supersonic self-sustained waves rely on 
the coupling between energy release in the media and the 
resulting mechanical waves (elastic, compression waves, 
etc.). The latter class of waves will henceforth be called 
detonations, in analogy to the self-sustained supersonic 
waves commonly encountered in reactive gases, combus- 
tion of dust particles in air and condensed phase ener- 
getic materials |T]. Such detonations can also appear 
in media sustaining thermo-nuclear reactions, and have 
been hypothesized to be the decomposition mode of white 
dwarfs undergoing supernova explosions [2^. Such super- 
sonic self-sustained waves have also been observed in a 
wide variety of elastic excitable media, such as Burridge- 
Knopoff models of frictional sliding, electronic transmis- 
sion lines, and active optical waveguides [3 . Likewise, 
phase change waves ^ , traffic jams [5 and shallow water 
waves [6^ also share the same characteristics of detona- 
tions, whereby the arrival of mechanical waves induces 
a change in the material state, conducive to the release 
of energy, which in turn modifies, or sustains the wave 
motion. 

A simple and elegant model for self-sustained super- 
sonic reaction waves has been introduced by Fickett [3|8]. 
This model is known to reproduce qualitatively many dy- 
namic traits of chemical detonations, such as the struc- 
ture of the self-sustained wave, initiation transients and 
response to boundary losses (see [8 ). Its mathematical 
simplicity offers a much simpler framework to study det- 
onations than the reactive Euler equations, as is the case 
for chemical detonations. Its simplicity also permits to 
consider this model as a unifying model to the wide va- 
riety of waves in excitable media mentioned above. The 
present study is an investigation of its spatio-temporal 



non-linear dynamics, which have not been addressed in 
previous studies. 

Indeed, in gaseous chemical detonations, experiments 
have demonstrated the large spatio-temporal variability 
in the wave propagation. In multiple dimensions, multi- 
scale cellular patterns have been observed [IJ, while in a 
single space dimension, a pulsating instability has been 
observed [9]. Numerical simulations of the reactive Eu- 
ler equations used to model one-dimensional pulsating 
detonations have shown the universality in the wave dy- 
namics [10 . As the sensitivity of the reaction rates is 
increased, stable travelling waves become oscillatory, and 
subsequently develop a hierarchy of period doubling bi- 
furcations appearing according to Feigenbaum's scaling 
pT] . and eventually become chaotic. At present, because 
of the complexity of the reactive Euler equations and 
resulting dynamics, neither the mechanism of the one- 
dimensional pulsating instability nor the reason for the 
universality in the period-doubling detonation dynamics 
are understood. In this paper we study the wave stability 
predicted by Fickett's model whose simplicity allows for 
analytical investigation of period-doubling and chaotic 
dynamics. 

Fickett's mathematical model for detonations is an ex- 
tension of the inviscid Burgers' equation to the reactive 
case, that is: 

dtp^d,p = (1) 
dtXr = r{p,Xr) (2) 

where all fields are defined on (x, t). The model is formu- 
lated in Lagrangian coordinates, where x denotes a mate- 
rial coordinate, while t represents time [8 . The variable p 
has the meaning of density in the reactive analogue. The 
flux term p in ([T]) has the meaning of pressure, see [8]. 
For simplicity, we choose the generic equation of state 
proposed by Fickett, that is 

P=l{p^ + KQ) (3) 
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The available energy in the excitable medium is Q and 
denotes the progress variable of the energy process, 
ranging from (un-reacted) to 1 (reacted). Equation ([2| 
describes the evolution of the energy release progress for 
each particle with Lagrangian coordinate x; a model for 
this reaction term will be described below. Note that by 
setting Q to zero, one recovers the well-studied inviscid 
Burgers' equation [12 . 

More insight into the interplay between wave motion 
and energy addition can be obtained by recognizing that 
the system of equations ([T]) and ^ is hyperbolic. It can 
be shown that the characteristic form can be written as: 
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From Q , we deduce that the system exhibits waves prop- 
agating forward with speed dx/dt = p. The wave commu- 
nicates changes in pressure amplitude only in the positive 
X direction. The amplitude of the wave varies with the 
heat addition Q at the rate r. Hence the model describes 
the physical property that waves may amplify in the pres- 
ence of energy release according to Rayleigh's criterion, 
i.e., if the energy release is in phase with the wave mo- 
tion. The second family of characteristics in ([5| gives the 
rate of energy release along a particle path. The physi- 
cal picture emerging is thus the reactivity set out along 
particle paths modifies the strength of waves propagat- 
ing forward. Through the coupling of the reaction rate 
(which we will ascribe below) to wave strengths, the feed- 
back loop is closed. Note that contrary to the reactive 
Euler equations used to model chemical detonations in 
fiuids, which admit three sets of waves [T3j, the analogue 
predicts two, as rear facing pressure waves are absent. 
This fundamental simplification permits to analyse the 
detonation problem in an analytically tractable system 
of equations, unlike the reactive Euler equations. 

The system admits a coherent self-propagating trav- 
elling wave solution having the properties of a detona- 
tion |8] . Although the details are available [8 , we briefiy 
describe its steady solution, as it serves as our start- 
ing point for the stability analysis. We seek a travel- 
ling wave solution to the system given by Q and 
The speed of the wave D can be found in terms of the 
unreacted state (po^^ro) in front of the wave and the 
reacted state (p2,A^2) behind the wave. For simplicity, 
and without any loss of generality, we set po = 0, A^o = 
and Xr2 = 1 to model an irreversible exothermic reaction. 
We also let p2 variable (i.e. the piston problem, see Fick- 
ett & Davis [1]). Adopting the notation [(] = C2 — Co, 
the resulting wave speed can be found (see [12 ) from 

D = mp] = iP2'+Q)/i2p2). 

The self-sustained travelling wave solution corresponds 
to the case where the forward propagating characteristic 
trailing the wave cannot penetrate the wave structure. 



and essentially represents an event-horizon. The speed 
of this so-called limiting characteristic thus needs to be 
equal to the detonation speed. Denoting this special case 
as the Chapman- Jouguet case (by analogy to the termi- 
nology used in chemical detonations [T) with subscript 
CJ, we require that p2 = D = Dcj^ from which we 
obtain the CJ speed of the detonation. 



Dcj = VQ 



(6) 



Since we are dealing with an inviscid system and the 
medium develops shocks according to Burgers' equation, 
the detonation can be assumed to be lead by an inert 
shock, across which there is no energy release and the 
density and pressure change discontinuously. We will de- 
note the state behind the shock with a subscript 1 (known 
as the von Neumann state in the chemical detonations). 
For a non-reactive shock satisfying the weak form of the 
inert inviscid Burgers equation, we obtain pi = 2D. 

The structure of the detonation wave across which en- 
ergy is deposited at a finite rate is obtained by integrat- 
ing the governing equations. The steady wave solution 
can be obtained by first introducing the change of co- 
ordinates = X — Dcjt — xo^t' = t) which defines a 
local coordinate system moving with the steady detona- 
tion. Making the formal change of variables and setting 
the time derivatives equal to zero in order to obtain the 
steady solution, we obtain: 



dC 
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XrQ] = 



(DcjXr) 



(7) 
(8) 



This system is integrated from the shock, with the inert 
shock state p = pi and A^ = as boundary condition at 
C = 0, once the rate r(p, A^) is given. 

In the present work, we propose and investigate a re- 
action model that is sufficiently simple to permit us to 
explain the unsteady period doubling wave dynamics of 
detonations and sufficiently rich to capture the main non- 
linearity in travelling waves in excitable media, which is 
the coupling between (shock) wave motion and exother- 
micity induced by the shock. We thus extend the models 
introduced in [8] and [14^, for which instability was pre- 
dicted from linear stability analysis, to a simple generic 
two step model. Following the leading shock, we as- 
sume the existence of a thermally neutral induction delay, 
whose duration depends on the strength of the shock. 
This is an excellent assumption for activated chemical 
reactions [1^, but can represent the excitability of any 
medium. Following the induction process, we assume an 
exothermic reaction that proceeds at a state-independent 
constant rate. The latter choice was selected in order to 
decouple the activation of the reactions with the dura- 
tion of the exothermic stage, in order to clearly isolate 
the physical phenomena governing wave motion and the 
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FIG. 1. Shock amplitude evolution for different values of a. 




instability mechanism and avoid the singularity in [8 as- 
sociated with square wave detonations. The resulting 
generic induction-reaction model we propose is thus: 



FIG. 2. Bifurcation diagram of the local maxima of detona- 
tion wave strength pmax for different induction zone sensitiv- 
ities a. 



dt\ = -K,H{\)e''y-^cj ^) (9) 

dt\r=Kr{l-H{\,)){l-\rY (10) 

where Ki and Kr are constants controlling the times 
scales of the induction and reaction zones, respectively. 
The Heaviside function i^(-) controls the timing of the 
onset of the second exothermic reaction, which starts 
when the induction variable reaches 0. Ahead of the 
shock, \i = 1 and = 0. We are also assuming that the 
reactions are only activated by the passage of the inert 
leading shock. The system to be solved is thus ([T]), ([9| 
and (ffol). 



The reaction model allows for direct analytical deriva- 
tion of the steady travelling wave solution. Ahead of 
the wave in the quiescent zone, we have, ^ > 0^ P = 0^ 
\i = 1 and A^ = 0. The induction zone terminates 
at Ci = —Dcj/Ki. In the induction zone we have 
Ci < C < 0, p = pi = 2Dcj, A, = 1 + K,/{DcjQ, 
and A^ = 0. For a reaction order v <1 the reaction layer 
terminates at a finite distance from the shock given by 
Cr = Ci ~ Dcj/{Kr{l — ly)). In the reaction layer, we 
have, p = Dcj (l + (1 + (1 - v)Kr/Dcj{C - Ct))^) 

and = 1 - (1 + (1 - v)Kr/Dcj{C, - Q)^ . 

The non-hnear stability of the travehing wave solution 
of 0, (|9| and ([lO) was investigated by numerical inte- 
gration starting with the steady travelling wave struc- 
ture as initial condition. The numerical integration uses 
the fractional steps method, whereby the hydrodynamic 
evolution and reactive step can be performed separately. 
The hydrodynamic step uses an exact first order Rie- 
mann solver. Owing to the simplicity of the reactive 
model, the reactive part of the governing equations is 
solved in closed form at each time step. A grid reso- 
lution of 256 grid points per detonation wave thickness 
was used, which ensured the stability boundary was grid 
independent with an accuracy of ±0.1 in the value of 
a. The results presented are for parameters, (5 = 5, 
i^i = 1, Kr = 2 and v = 0.5. Below the critical value 



a = 5.7, the steady solution is stable, and propagated 
with the steady wave structure given above at its con- 
stant CJ speed given by (|6|. Above this critical value 
the travelling wave solution is unstable and develops a 
stable limit-cycle, as shown in Fig. [l] As a increases, the 
amplitude of the pulsations also increases until a period 
doubling bifurcation occurs at a = 6.9. Fig. |2] shows 
the bifurcation diagram. Further increases in a yields 
another bifurcation at a = 7.7, followed by subsequent 
bifurcations occurring with smaller changes in a. Eventu- 
ally, the pulsations become chaotic, with isolated values 
of a for which the dynamics have an odd period; one ex- 
ample is shown in Fig. [l] This implies the onset of chaos 
[10 . The sequence of period doubling bifurcations and 
onset of chaos is very similar to the non-linear dynamics 
of chemical detonations |T0] and many other non-linear 
systems. The results thus clearly highlight that our sim- 
ple detonation model captures this universality observed 
in non-linear dynamics of complex systems. 

The analogue model also allows for a straightforward 
interpretation of the instability mechanism in detona- 
tions. In order to study the non-linear instability mech- 
anism of the proposed detonation analogue, we focused 
our attention on the single mode limit-cycle obtained for 
a = 6.8. Figure |3] illustrates the evolution of the wave 
structure over approximately two oscillation periods, in 
the frame of reference of the steady travelling wave. To 
visualize the dynamics, we reconstructed an (arbitrary) 
discrete set of pressure waves by integrating Q starting 
from arbitrary locations. We used a predictor-corrector 
method and interpolated on the solution obtained above. 
The lead shock front of the detonation corresponds to 
the locus where these characteristics coalesce. Behind 
the oscillating lead shock are the two zones of induction 
and reaction. Note that the finite (numerical) dissipation 
at the discontinuity at the rear of the reaction zone makes 
the characteristics bend somewhat towards the reaction 
zone in a very narrow region. Away from this region. 
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FIG. 3. Space time diagram illustrating the pressure waves 
in the reaction zone of a pulsating detonation. 



by virtue of the characteristic equation Q, the pressure 
waves have constant amplitude and speed everywhere ex- 
cept in the reaction zone, where they accelerate owing 
to the heat release. By investigation of the character- 
istic diagram of Fig. |3) the detonation wave structure 
can be easily interpreted as the coherent wave structure 
formed by the amplification of forward travelling waves. 
These are amplified across the reaction zone and even- 
tually reach the hydrodynamic shock. Since the onset of 
the reactions is controlled by the lead shock, the pres- 
sure waves continuously see the same reacting field and 
the self-sustained detonation phenomenon occurs. 

The instability mechanism itself can be inferred from 
the characteristic diagram shown in Figure [3j First, one 
can note that the oscillations of the leading shock pro- 
vide a modulation in the duration of the induction zone 
and location of the reaction zone. This is due to the in- 
duction time sensitivity on shock state given in ([9|. As 
can be verified in Fig. |3j the points when the shock is 
strongest (S) correspond to the shortest induction times, 
which is communicated along the particle paths (P). Like- 
wise, the weakest shocks (W) yield the longest induction 
times. Due to this modulation in the onset of exother- 
micity, the acceleration and deceleration phases of the 
lead shock pulsations can be deduced. The amplifica- 
tion stage of the lead shock corresponds to the arrival 
of pressure waves that travelled in phase with the energy 
release zone, which can be identified as the regions where 
the pressure waves travel almost parallel to the reaction 
zone band. Likewise, the deceleration phases correspond 
to the arrival at the lead shock of waves travelling out 
of phase with the energy release zone. The waves travel- 
ling in phase with the energy release amplify more, and 
communicate an acceleration to the lead shock. Since 
this occurs during the lead shock amplification stage, the 
feedback accentuates the amplification. Likewise, a de- 
celerating shock provides a non-coherent interaction be- 



tween the forward pressure waves and exothermicity, fur- 
ther promoting the deceleration. In our system, the co- 
herent amplification can be obtained by integrating Q 
for a constant rate, which shows that the amplification of 
a pressure wave is proportional with the residence time 
of the wave in the energy release zone. When a pressure 
wave stays in phase with an energy release zone for longer 
times, it acquires the most amplification. The pulsation 
mechanism, which controls the sequence of acceleration 
and deceleration phases can also be seen in Fig. [Sj Fol- 
lowing an acceleration stage, the forward characteristics 
emanating from the rear of the reaction zone, which only 
clip the reaction zone and obtain very little amplification 
form the expansion waves (E). The expansion waves im- 
mediately following the compression waves provide the 
restoring mechanism for the instability. Note that the 
same mechanism has been suggested to be at play in 
chemical detonations [13[ [15], although the complexity 
of the governing equations did not permit to clearly iso- 
late these effects. The much simpler detonation model 
suggested in the present study provides very similar dy- 
namics and permits a much more accurate physical in- 
vestigation of the physical mechanisms controlling the 
instability, namely the coherent amplification of forward 
facing waves modulated by the onset of reactions. This 
is essentially Rayleigh's criterion for (acoustic) wave am- 
plification by coherent energy release. 
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